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We present a method to calculate the electronic charge density of periodic solids in the GW approx- 
imation, using the space-time method. We investigate for the examples of silicon and germanium 
to what extent the GW approximation is charge-conserving and how the charge density compares 
with experimental values. We find that the GW charge density is close to experiment and charge 
is practically conserved. We also discuss how using a Hartree potential consistent with the level of 
approximation affects the quasi-particle energies and find that the common simplification of using 
the LDA Hartree potential is a very well justified. 



I. INTRODUCTION 

Amongst the established methods to calculate the elec- 
tronic properties of solids, Hedin's GW approximation!!] is 
notable for its unrivalled success in predicting the energy 
gaps of semiconductors and insulators (e.g. Refs. ||-|5|). 
Unlike alternative methods which are based on density 
functional theory and describe ground state properties 
by mapping the many-electron problem onto an effective 
one-electron problem, it seeks to find a solution for the 
exact one-particle Green's function of the many-electron 
system to which it is applied. With this knowledge, a 
wide range of properties of the system under considera- 
tion can be accurately described, among them the opti- 
cal excitation spectrum, the density of states, the charge 
density and the total energy. 

As an exact solution of the many-electron problem re- 
mains today as elusive as ever, any attempt to determine 
the Green's function has to rely on skilled approxima- 
tions. The problem of finding the exact Green's function 
of a system is equivalent to finding its self-energy, fn the 
GW approximation, this quantity is approximated by the 
product of G, the Green's function, and W, the screened 
Coulomb interaction; hence the name. 

GW calculations for semiconductors so far have con- 
centrated on the optical excitation spectrum, the one 
area where GW has had its most striking successes, but 
other aspects of electronic structure have not been in- 
vestigated with this method for real materials. This is 
largely due to the prohibitively large numerical effort. 
Recently, however, a new technique has been developed, 
the GW space-time method,til that makes GW calcula- 
tions much faster and allows a precise computation of the 
energy dependence of the self-energy without recourse to 
plasmon pole models. This opens up the way to new 
applications for the GW formalism. 

In this paper we will concentrate on the electronic 
charge density^ of semiconductors, present a new tech- 
nique for extracting it from the Green's function at GW 
level and show for the first time results of GW charge 
density calculations for real materials at the example of 
Si and Ge, for which very careful analyses of the experi- 



mental data are available for comparison .Effl We also in- 
vestigate the problem of charge conservation and the in- 
fluence of using a Hartree potential consistent with the 
level of approximation on the quasi-particle spectrum. 

Atomic units are used throughout unless otherwise 
noted. 



II. THEORY 

The central equation of the GW approximation is the 
expression for the self-energy: 

E(r, r'; r) = iG(r, r'; r)W(r, r'; r). (2.1) 

Full self-consistency in the GW approximation would 
mean that the Green's function that enters into Eq. (2.1) 
were itself a solution of the equation 



-^V 2 + V ext (r) + V£ c {r) + E 5c (r, r») G(r, r', u) 

= -iS(r,r') (2.2) 

where the superscript SC for S and the Hartree potential 
Vh indicates that these quantities have been determined 
self-consistently. In the case of the Hartree potential this 
means the charge density at GW level has to be known. 
Such a self-consistent GW calculation has nevee. been 
done for a real system and there are indicationalj that 
self-consistency will destroy the good agreement between 
the computed and experimental quasi-particle spectrum. 

In this paper we employ the usual one-iteration non- 
sclf-consistent GW approximation. However, we will de- 
termine the charge density and thus the Hartree potential 
self-consistently at a level where the self-energy is kept 
fixed at its first-iteration value. The charge density itself 
is an indicator as to whether such an approach is jpsjt-i-i 
fied, as it has been shown by Baym and KadanofOE 2 ] 
that generally only self-consistent approximations to the 
self-energy conserve particle number strictly. We will 
therefore monitor carefully whether our approach vio- 
lates particle number. 
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The traditional way to set up and solve the GW equa- 
tions has been to express and compute all quantities in 
the reciprocal space and energy domain. This involves 
sums scaling with the fourth power of the number of 
plane waves N used tcurepresent the wavefunctions. Ro- 
jas, Godby and NeedsQ have shown recently that using 
the most appropriate representation at the various stages 
of the computation, either real space and time or recip- 
rocal space and energy and changing between represen- 
tations by means of Fast Fourier Transforms, can bring 
this scaling down to iV 2 , thus leading to very significant 
time savings. 

Another point noted by Rojas et al. is that a detailed 
calculation of all relevant quantities along the imaginary 
time or imaginary energy axis is possible without recourse 
to a plasmon pole or similar model. This is because on 
the imaginary time axis the Green's function, the po- 
larisability and the self-energy are all rapidly decaying 
smooth functions which can easily be represented on a 
regular grid. Transformation to imaginary energy via an 
FFT is straightforward. If the properties of £ on the 
real energy axis are needed, analytic continuation with 
the help of model functions is possible. However, for the 
calculation of the charge density alone, the knowledge 
of the self-energy on the imaginary energy axis only is 
sufficient. 

Having determined the self-energy £ on the imaginary 
energy axis as described in Ref. ^, with 



E(r, r'; it) = *G (r, r'; %t)W(t, r'; *r), 



(2.3) 



where Go is the Green's function at LDA level, the 
Green's function G at GW level obeys the Dyson equa- 
tion 



G(iuj) — Gq(ilu) + Go(iw) x 

x (E(iw) + AV H - V xc - e ) G(iw) 



(2.4) 



where G, Go, £, AVh and V xc in this notation are to 
be understood as matrices in a plane-wave basis and ma- 
trix multiplication of the factors on the right hand side 
is implied. AVh is the change in the Hartree potential 
due to the density change and has to be determined self- 
consistently, V xc is the LDA exchange correlation poten- 
tial and eo represents the shift in the Fermi level with re- 
spect to the vacuum in the GW calculation from its value 
in theADA. The use of this shift eo was first proposed by 
Hedintl but is often neglected in GW calculations. We 
use it here to introduce an element of self-consistency 
into the equations by ensuring that the Fermi energy is 
the s ame before and after applying the GW correction. 
Eq. (2.4) can directly be solved by matrix inversion. 

The relationship between charge density p and Green's 
function is given by the equation 



2 f° 

* J-oc 



(2.5) 



Since we know the charge density at the LDA level of 
approximation we only need to evaluate the charge den- 
sity difference Ap between the LDA and the GW result. 
The Green's function has poles in the second and fourth 
quadrant of the imaginary plane, just infinitesimally off 
the real axis. Using Cauchy's Theorem and the fact that 
for complex energy z AG(z) vanishes quadratically as 
\z\ — > oo, where AG = G — Go, we know that the charge 
density difference can just as well be calculated by an 
integration along the imaginary energy axis: 



2 r° 

Ap(r) = / dujRcAG(r,r;iuj). 

7T J-oo 



(2.6) 



This means that the knowledge of the Green's func- 
tion along the imaginary axis is sufficient to calculate 
the charge density at GW level without the need for ex- 
plicit analytic continuation to the real energy axis, which 
would inevitably have to involve some sort of model. 

In theory, the integration of AG(iuj) according to Eq. 
(2.6) is straightforward. In practice, in order to achieve 
satisfactory convergence in a numerical integration, one 
would have to include such a large number of energy 
points in the integration, as to make it prohibitively ex- 
pensive, because for each energy point the self-energy 
would first have to be computed and the Dyson equation 
solved. On the other hand, we know that the self-energy 
as a function of imaginary energy decays quadratically 
for large energies. This allows us for large u> to make an 
expansion of G(iu>) in powers of l/u> and neglect terms 
of order higher than 1/lu 2 . The perturbatively treated 
high-energy tail can be integrated analytically from some 
energy loq to infinity, as we wi ll show in the appendix. 
Up to luq we integrate Eq. ( |2.6| ) numerically, loq is thus 
treated as a convergence parameter. 

For eac h new An we compute a new AVh and solve 
Eqs. (2.4) and fl2.6| ) repeatedly until AVh is stable. 



where it is assumed that the zero of energy is chosen at 
the Fermi level. 



III. RESULTS 

Before we can make comparisons to experimental re- 
sults we have to account for the fact that we work in 
the pseudopotential approximation. This means that the 
charge density computed as described in the previous sec- 
tion contains no contribution from the core electrons and 
the contribution of the valence electrons is modified in the 
core region as well. ■— ■ 

We follow Nielsen and MartinCJ in writing the total 
charge density as the superposition of the atomic densi- 
ties plus a deformation density pdef, which is defined in 
reciprocal space as 

Pdof(G) = p S olid : ps(G) — S(G)p a tom,ps(G), (3.1) 

where G is a reciprocal lattice vector, S the structure 
factor of the crystal and /9 so iid,ps the charge density gen- 
erated by the pseudowavefunctions of the valence elec- 
trons in the solid and p a tom,ps the same in the free atom. 
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The Fourier coefficients of the total charge density in the 
solid are then computed according to 



Asolid(G) = 5(G)p atom (G) + Pdef(G) 



(3.2) 



where p a tom is the full charge density of valence and core 
electrons of a free atom. This approach assumes rigid 
cores, in line with the assumptions underlying the pseu- 
dopotential method. 

Where we find that the charge in the GW approxima- 
tion is not strictly conserved, we normalise the density 
by multiplying all structure factors computed according 
to Eq. (3.2) by a constant that assures the overall correct 
number of electrons per unit cell. 

The results we present in this section come from LDA 
calculations that were performed with pseudopotentials 
generated by the Hamann method,La a plane-wave cut- 
off of 17Ry for Si and 20Ry for Ge and 10 special k 
points. For Ge relativistic effects were included in the 
pseudopotential but no spin-orbit coupling was taken into 
account in the solid. A Ceperley-Alder exchange cor- 
relation potentiallia rwas used in the parametrisation by 
Perdew and ZungerJl3 In the GW space-time calculation 
we used 65 bands truncated to 169 plane waves from 
the LDA calculation. The time grid comprised 60 points 
spaced at 0.314 a.u. which were zero-padded to 120 points 
before transformation to imaginary energy, giving an loq 
of 10 Hartrees. Other parameters of the GW calculation 
were as described in Ref. ||. 

Table || shows our computed structure factors of the 
pseudo valence charge of silicon for several reciprocal 
lattice vectors G. The second column shows the LDA 
values and the third the GW correction Ap. Ap is ac- 
tually the difference between the normalised GW den- 
sity and the LDA density. Without normalisation the 
GW valence density integrates to a total of 8.0233 elec- 
trons per unit cell (instead of 8), a charge violation of not 
quite 0.3%. This observation of a very small but finite 
charge violation is in line with rigocous analytical results 
by Schindlmayr for a model systemO The fourth column 
lists the (normalised) GW correction to the density that 
results if the Hartree potential is kept fixed at its LDA 
value. As one would expect, this unscreened correction is 
larger. However, because contributions from all energies 
are integrated over and we are looking at non-zero G the 
screening effect is much weaker than one might naively 
expect by scaling down with the dielectric constant. The 
LDA structure factors of the pseudo valence charge are 
themselves not very meaningful, as they depend to some 
extent on the specific pseudopotential used. They are 
only listed to give an idea of the relative magnitude of 
the GW correction. 

Table || shows the structure factors of the total den- 
sity of Si for several reciprocal lattice vectors. The sec- 
ond column lists the LDA values and the third the (nor- 
malised) GW values. The biggest contribution comes 
actually from the cores, as the comparison with the va- 
lence structure factors in Table | shows, so that now the 



difference between the LDA and the GW structure fac- 
tors looks even less significant. Because of the weight 
of the core contribution the lattice constant has a cru- 
cial influence on the structure factors. We have used for 
both the LDA and the GW calculation the experimental 
value of 5.43A. In the last column we list experimental 
static (zero temperature) structure factors for compiarA 
son. They are the Fourier transforms of a model fitlij' 
to experimental data. The agreement of both GW and 
LDA with experiment is very good. 



Table III gives the valence structure factors of germa- 
nium and the GW correction to it. The raw GW valence 
density integrates to 8.004 electrons per unit cell (instead 
of 8), a charge violation of 0.05%. The values shown are 
for the normalised density and values for the unscreened 
correction are given as in Si. Table |v| compares the full 
LDA and GW density structure factors in bulk Ge with 
experimental static structure factors from Ref. |^ (fit C 
in their Table III). The agreement of both LDA and GW 
values with experiment is again very good. We have used 
a Ge lattice constant of 5. 66 A in LDA and GW calcula- 
tions. 

Note that the structure factors listed in the tables have 
to be multiplied by the diamond lattice structure factors 
5(G) to get the Fourier coefficients of the charge density 
in the crystal as expressed by Eq. (3.2). 

Small though the differences between the LDA and 
GW densities are, we can see that for both Si and Ge 
the GW charge in real space is slightly less concentrated 
near the atom sites and moves a little into the intersti- 
tial and bonding regions. As a trend this goes into the 
right direction since it is known that the LDA tends to 
accumulate too much charge near the atoms. Figs, [l] and 
H show the charge density differences in Si and Ge along 
the edge of a conventional cell between two atom sites. 

We give a list of quasiparticle energies at several sym- 
metry points for Si in Table [v] and for Ge in Table VI . 
The second column gives the values with the Hartree po- 
tential taken at LDA level and the third with a Hartree 
potential which uses the GW density. One can see that 
the influence of AVh on the quasi-particle energies is very 
small. 



IV. CONCLUSION 

The results presented in the previous section show that 
in the materials we have investigated the GW approxi- 
mation can for most practical purposes be considered to 
be charge conserving. It could be demonstrated that the 
charge density is very close to experimental values and 
also little different from the density at LDA level. This 
gives support to the common practice of computing the 
quasi-particle corrections to the LDA eigenvalues without 
adjusting the Hartree potential. Adjusting the Hartree 
potential to the GW density has only a small effect on 
the quasi-particle energies. 
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V. APPENDIX 

To show ho w w e treat the high-energy tail of the inte- 
grand in Eq. (2.6) let us rewrite Eq. ( |2.4| ) in the form 



G(iu>) = F-\iuj)Go(iuj), 



(5.1) 



F(iw) = (1 - GoM [E xc (»w) + AV H - V xc - e }). 

(5.2) 

For large u, as Go decays as 1/w, we can expand F^ 1 
in powers of Go, keeping in mind that £ itself contains 
one factor Go- To first order in Go then 

lim F-\ioj) = 1 + GoM + Ay H - V xc - eo] , 

U) — > — oo 

(5.3) 

where T, x is the exchange part of the self-energy that does 
not itself depend on the high energy tail of Go, 



EsM = G o (0)V c , 



(5.4) 



where V c denotes the inter-electronic Coulomb interac- 
tion. 

To this order of approximation then 
lim AGUlj) 

UJ— i-OO 

= G (iu>)lZ x + AV H -V xc \G (iuj). (5.5) 
After a few transformations we can find from this 
lim ReAG(r, r; ioS) — 

LJ — '' — OO 

E { RC [*nk(r)*n'k(r)Pnn' (k)] /$(<«/; V 



k Tin' 



-Re 



[*nk(r)*„' k (r)P nn ,(k)] /S(«;k)} (5.6) 



where 



/$( W ;k) = Re ' 



(w nk - «w)(w„/ k - ' 



/f,(u,;k)=Im ' 



(Wnk - iu>)(Un'V. ~ U*)) ' 



(5.7) 



(5.8) 



P w (k) = E *„k(G)[£ x (k,G,G') 

G,G' 

+ AVfc (G - G') - V XC (G - G')] *„'k(G'). (5.9) 

The advantage of this formulation is that the en- 
ergy dependent parts can be integrated an alyt ically. We 
therefore split up the integration in Eq. (2.6) into two 
parts, the first of which is performed numerically up to 
some suitably chosen value loq and the remainder of the 



integral from u>o to — oo is evaluated analytically by using 
the relationship 



/«(k)H^ 



_ j "i,k-"„k' 

| 1 u) Q /u) n 



arctan — arctan -^ SL - 



OJ„k l + (wo/"„k) 2 



if W nk 7^ w "'k 
if w nk = w n'k 



and 



^7^T ln ifWnk^^n'k 



if w„k = w ra 'k 



[ _ 2(t l ;> 2 +t l ;< 2 ) 

w< = min(|w nk |, |w nk '|), w> = max(|w nk |, |w n k' 
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TABLE I. Structure factors (absolute values) of the 
pseudo valence charge density of Si [e/atom] at LDA level and 
the correction from the GW approximation with and without 
adjusting the Hartree potential. 



G 



LDA valence 



Ap 



Ap 



4 







(V H GW) 


(V„ LDA) 


1 1 1 


1.2487 


-0.0054 


-0.0076 


2 2 


0.0323 


-0.0035 


-0.0043 


Q 1 1 

oil 


n o a on 




n nnoQ 


4 


0.1877 


0.0009 


0.0012 


3 3 1 


0.0601 


0.0021 


0.0023 


4 2 2 


0.0685 


0.0008 


0.0010 


3 3 3 


0.0755 


-0.0003 


-0.0001 


TABLE II. Structure factors (absolute values) of the total 


charge 


density of Si [e/atom]. 






G 


LDA 


GW 


Exp. 


1 1 1 


10.7210 


10.7157 


10.713 


2 2 


8.6536 


8.6501 


8.655 


3 1 1 


8.0205 


8.0182 


8.027 


4 


7.4414 


7.4423 


7.454 


3 3 1 


7.2256 


7.2277 


7.246 


4 2 2 


6.6984 


6.6992 


6.712 


3 3 3 


6.4086 


6.4083 


6.420 


TABLE III. Structure factors (absolute 


values) of the 


pseudo 


valence charge density of Ge [e/atom] at LDA level 


and the correction from the GW approximation with and 


without adjusting the Hartree potential. 




G 


LDA valence 


Ap 


Ap 






{V H GW) 


(V H LDA) 


1 1 1 


1.3181 


-0.0169 


-0.0239 


2 2 


0.0035 


-0.0086 


-0.0107 


3 1 1 


0.2201 


-0.0036 


-0.0033 


4 


0.2052 


0.0025 


0.0037 


3 3 1 


0.0974 


0.0037 


0.0040 


4 2 2 


0.0865 


0.0027 


0.0033 


3 3 3 


0.0726 


0.0017 


0.0023 





TABLE V. 


Quasi-particle 


energies of Si [eV]. 




level 


T T~~» \ 

LDA 


GW + V H 


GW + V H 


Exp. 


Tic 

^25'v 

ri5 C 
rv c 


-11.88 
0.00 
2.59 
3.26 


-11.91 
0.00 
3.26 
4.05 


-11.91 
0.00 
3.26 
4.03 


-12.50 
0.00 
3.05 
4.1 


Xl v 
Xav 
X\ c 
Xac 


-7.77 
-2.81 
0.62 
10.11 


-7.88 
-2.92 
1.24 
11.00 


-7.88 
-2.92 
1.25 
10.99 


-8.18 
-2.9 
1.25 
10.95 


L^v 
Liv 
Ly-u 
Lie 
Lzc 
Ly c 


-9.56 
-6.95 
-1.16 
1.46 
3.34 
7.73 


-9.64 
-7.09 
-1.22 
2.14 
4.07 
8.34 


-9.64 
-7.08 
-1.22 
2.14 
4.08 
8.36 


-9.3 
-6.7 
-1.2 
1.65 
4.15 


Gap 


0.49 


1.10 


1.11 


1.17 




TABLE VI. 


Quasi-particle energies of Ge [eV]. 




level 


LDA 


GW + V L H DA 


GW + V% w 


Exp. 


Tic 

rv c 


-12.63 
0.00 
0.01 
2.61 


-12.79 
0.00 
0.66 
3.11 


-12.79 
0.00 
0.63 
3.14 


-12.60 
0.00 
0.89 
3.21 


Xl v 
Xav 
Xlc 
Xzc 


-8.83 
-2.98 
0.68 
9.56 


-8.97 
-3.11 
1.08 
10.32 


-8.98 
-3.09 
1.12 
10.32 


-9.3 
-3.15 
1.3 


L^v 
Liv 
L-i'v 
Lie 
Lzc 

L2'c 


-10.61 
-7.54 
-1.35 
0.14 
3.77 
7.24 


-10.76 
-7.69 
-1.40 
0.64 
4.29 
7.65 


-10.76 
-7.68 
-1.40 
0.66 
4.32 
7.73 


-10.6 
-7.7 
-1.4 
0.74 
4.3 
7.8 



TABLE IV. Structure factors (absolute values) of the total 
charge density of Ge [e/atom]. 



G LDA GW Exp. 

1 1 1 27.5205 27.5036 27.453 

2 2 23.6819 23.6734 23.677 

3 1 1 22.1675 22.1639 22.138 

4 20.3205 20.3230 20.273 

3 3 1 19.4545 19.4582 19.509 

4 2 2 18.0361 18.0388 18.066 
3 3 3 17.2916 17.2939 17.315 



FIG. h Difference between the GW and LDA charge 
densities^ of Si (atomic units) along the edge of a conventional 
cell between two atom sites. Distance in Bohr radii. The GW 
correction slightly increases the concentration of electrons be- 
tween the atoms. 

FIG. 2i Difference between the GW and LDA charge 
densitiesQ of Ge (atomic units) along the edge of a conven- 
tional cell between two atom sites. Distance in Bohr radii. 
The GW correction slightly increases the concentration of 
electrons between the atoms. 



5 



0.0004 
0.0003 
0.0002 
0.0001 


-0.0001 
-0.0002 



0.001 




~0 0004 1 1 1 ' 1 1 1 ' 1 1 1 ' 1 1 1 ' 1 1 1 ' 

2 4 6 8 10 

x (a.u.) 



